AD-A102  943 
UNCLASSIFIED 


OHIO  STATE  UNI V  COLUMBUS  DEPT  OP  GEODETIC  SCIENCE  F/®  8/5 

FEASIBILITY  STUDIES  FOR  THE  PREDICTION  OF  THE  GRAVITY  DISTURBAN— ETC (UJ 
MAR  81  H  SUENKEL  FI9628-79-C-0075 

DGS-311  AFGL-TR-81-0084  NL 


AFGL-Tr-81-0084 


LEVE# 


AFGL-TR-81-0084 


CO 


FEASIBILITY  STUDIES  FOR  THE  PREDICTION  OF  THE 
GRAVITY  DISTURBANCE  VECTOR  IN  HIGH  ALTITUDES 


Hans  SSnkel 


Department  of  Geodetic  Science 
The  Ohio  State  University 
Research  Foundation 
Columbus,  Ohio  43210 


March  1981 


DT!C 

r 


AUG  1  7iS81  j 


H 


Scientific  Report  No.  5 


Approved  for  public  release;  distribution  unlimited 


AIR  FORCE  GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 
HANSCOM  AFB,  MASSACHUSETTS  01731 


81  8  1?  05 


Qualified  requestors  may  obtain  additional  copies  from  the 
Defense  Technical  Information  Center.  All  others  should 
apply  to  the  National  Technical  Information  Service. 


Unci  ass  if ied 


SECURITY  CLASSIFICATION  of  this  PAGE  fWhon  Data  Enltrmd) 


1  fj  •  report  documentation  page 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

Iv-AVpoHT  NUI>»»-  .  Is  GOVT  ACCESSION  NO. 

i  AFGiy/TR. 81-0084  ]rfD-A/OH 

3  RECIPIENT'S  CATALOG  NUMBER 

9V3 

4.  TITLE  (and  Subtitlm) 

> feasibility  studies  for  the  prediction  of 

-  THE. GRAVITY  DISTURBANCE  VECTOR  IN  HIGH 
ALTITUOES^ 

5.  Type  op  report  a  perioo  covered 

Scientific  Report  No. 5 

6.  PERFORMING  ORG.  REPORT  NUMBER 

Report  No . 311 

*r  AOTHOaCa».  — -  m 

\  Hans/ Sl^ikel  / 

»  CONTRACT  OR  GRANT  NUMBER/*/ 

F19628-79-C-0075 

» 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Department  of  Geodetic  Science 

The  Ohio  State  University 

Columbus,  Ohio  43210  ^ 

MM 

II.  CONTROLLING  OFFICE  NAME  AND  AOORESS 

Air  Force  Geophysics  Laboratory  // 

Hanscom  AFB,  Massachusetts  01731  — ' 
Project  Monitor  -  Bela  Szabo/LW 

12.  REPORT  DATE 

March  1981 

I*.  WUMaew  0~F"^~AGES 

53  pages 

rt  MONITORING  AGENCY  name  4  AOORESS (II  dl/Uranl  Iron,  Controlling  Olllco) 

I'i 

'> 

15.  SECURITY  CLASS,  (ol  thi*  report) 

Unclassified 

IS*.  DECLASSIFICATION/  DOWNGRADING 

schedule 

16.  distribution  statement  c of  thi,  noport) 

Approved  for  public  release;  distribution 

uni imi ted 

~  „  nri  n 

17.  DISTRIBUTION,  STATEMENT  (of  thm  abstract  mntarad  In  Block  20,  It  diifarant  from  Report)  ;*  %  v  \  v. 

•7:"* 

.  __  A,  -  ; 

18.  SUPPLEMENTARY  NOTES 

*  T 

18.  KEY  WOPOS  (Continum  on  rararaa  aid*  it  n«c«i««ry  and  identity  by  block  numbmrj 

disturbing  potential 
least-squares  collocation 
integral  formula 
isotropic  integral  kernel 
frequency  domain 

20-  ABSTRACT  (Confirm*  on  r*v*r# •  aid*  It  nacaaaary  and  identity  by  block  tivmbmr) 

/DnTJnnacCo^cL2frthe,9ravit^  disturbance  vector  in  high  altitude 
(30.000  -  200-000  feet),  predicted  from  a  surface-cover i ng  set  of 
mean  gravity  anomalies,  is  estimated.  Two  methods  are  used  and 
found  to  provide  estimates  which  differ  by  less  than  10%,  the 
least-squares  collocation  and  the  integral  solution®  for  the  inte¬ 
gral  solution,  the  estimation  of  the  representation  error  has  been 
performed  in  the  frequency  domain.  For  the  collocation  solution 
an  optimal  algorithm  has  been  develoDod  whirh  tai-oc  ,i«<ujii  +  ,no  • 

00  1  JAH  TJ  1473  .-SOITION  OF  I  NOV  «S  IS  OBSOLETE  U  0  C  1  3  S  S  i  f  i  e  d 


✓  / 


■/ 


SECURITY  CLASSIFICATION  OF  THIS  PACE  fPhan  Dmlm 


SECURITY  CkASSinCATIOM  OT  This  RAGEtWiwi  Dmtt  Sam r*d) 


the  regular  data  distribution  and  is  up  to  64  times  faster  than 
a  non-optimized  solution.  "  '  ~ 

The  results  indicate  that  the  radial  component  of  the  gravity 
disturbance  vector  can  be  estimated  with  an  accuracy  of  ±1  mgal 
at  an  altitude  of  about  50-000  ft  on  the  basis  of  the  available 
data  sets;  in  order  to  achieve  the  same  accuracy  at  30-000  ft 
altitude,  the  data  error,  particularly  chat  of  5'  x  5‘  anomalies 
has  to  ,be  reduced  by  some  6055;  the  available  data  distributions 
are  adequate.  The  prediction  error  drops  quickly  with  increasin; 
al  ti  tude. 
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component:  with  the  best  available  data  distribution  an 
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(this  corresponds  to  ±0115  in  the  direction  of  the  gravity 
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a  spherical  distance  of  about  30°;  in  addition,  the  overall  data 
error  needs  to  be  reduced  by  some  30%.  The  prediction  error 
decreases  only  slowly  with  increasing  altitude. 
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1 .  PROBLEM  FORMULATION 

Given  a  set  of  mean  free  air  gravity  anomalies  at  sea  level 
and  associated  ms-error  estimates ,  the  accuracy  of  a  predicted 
gravity  disturbance  vector  in  high  altitude  (30  000  -  200  000  ft) 
should  be  estimated.  The  underlying  prediction  concept  is  least- 
squares  collocation  based  on  a  homogeneous  and  isotropic  global 
covariance  function  for  free  air  gravity  anomalies,  which  is  sup¬ 
posed  to  provide  best  linear  unbiased  estimates  with  "best"  depen¬ 
ding  on  the  choice  of  the  covariance  function  and  on  the  data  error 
characteristic . 

Although  a  complete  and  consistent  prediction  algorithm 
(a  sub-module  of  GSPP)  was  available  (Siinkel,  1980)  ,  the  design 
of  a  strictly  problem-oriented  procedure  turned  our  to  be  necessary: 
the  prediction  algorithm  of  GSPP  is  designed  for  the  general  case 
of  a  heterogeneous  set  of  irregularly  distributed  data;  it  does 
not  take  advantage  of  symmetries  in  the  data  distribution.  Since 
the  data  sets  to  be  studied  are  fairly  large  (  >  2000) ,  a  straight¬ 
forward  least-squares  collocation  solution  turned  out  to  be  prohibit¬ 
ive  because  of  excessive  computation  time  requirements.  Therefore, 
studies  have  been  carried  out  which  led  to  the  realization  of  a 
computer  routine, well  tailored  for  the  problem  in  consideration: 
blockwise  homogeneous  data  which  are  symmetrically  distributed  with 
respect  to  the  computation  point.  Depending  on  its  latitude,  the 
calculation  speed  can  be  increased  by  a  factor  of  at  least  8  and 
at  most  64  relative  to  non-problem-oriented  algorithms.  Tests  and 
comparisons  with  the  existing  prediction  module  have  been  carried 
out  with  small  data  sets. 

Stimulated  by  the  result  of  practical  experiments  (Rapp  and 
Agajelu,  1975)  showing  that  collocation  solutions  differ  from  Pois¬ 
son  integral  solutions  by  about  10  %  only,  detailed  studies  were 
performed  using  this  approach.  Rapp's  findings  are  largely  confirmed 
by  my  results;  this  is  quite  remarkable  since  the  integral  estimation 
presented  here  is  almost  exclusively  performed  in  the  frequency 
domain. 


THE  DATA 


The  data  set  to  be  used  in  these  studies  consists  of  mean 
free  air  gravity  anomalies  covering  the  surface  of  the  sphere 
partly  as  well  as  totally.  The  arrangement  of  the  data  is  symmet 
rical  with  respect  to  the  computation  point:  Starting  with 
5'  x  5'  mean  anomalies  in  a  rectangular  region  with  the  computa 
tion  point  P  as  its  center,  rectangular  zones  are  defined,  eac 
zone  being  completely  covered  by  mean  anomalies,  15'  x  15'  , 

1°  x  1°  and  5°  x  5°  means  (Fig.  2.1) . 


FIG.  2.1  Data  arrangement 


The  prediction  (computation)  point  is  assumed  to  have  a  latitude  of 
*»  40°  . 


As  far  as  the  data  distribution  and  the  a  priori  estimates 
of  corresponding  rms-estimates  are  concerned,  18  cases  had  to 
be  studied.  The  following  two  tables  summarize  these  cases. 
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Case 

fi  l  ( 5  '  x  5  ' ) 

*2  (15* x 15') 

a3(rxi°) 

/  r-  0  p  0  i  1 

^  x5  )  1 

af  ( 5  J  <  5  3 ) 

1 

2"  *2° 

6°  x  8° 

26  °  x  30° 

50° ' 70°  . 

180  "  < 360 

2 

3°  x  4  0 

-0  no 

7  x9 

26  0  x  30 

50°  x  70  3  i 

180''  •  360 ' 

3 

-  0  .0 

3  x4 

-0  ~  0 

7  x9 

30°  x  34° 

50° *70°  ; 

180°  x 360' 

4 

3°  x  4° 

7°x9° 

383  x  42° 

50°  x  70°  | 

1 803  x  360 

5 

*%  0  i  0 

3  x  4 

10°  x 1 2° 

26  0  x  30° 

50°  x  70°  ; 

1 

180°  <360° 

6 

-j  0  ,,  0 

3  x  4 

16°  x 1 8° 

26°x30° 

50°  x  70°  * 

1 80“ x  360 ' 

7 

7°x9° 

10°x 1 2° 

26  0  x  30° 

50c  x 70°  1 

1 80“ x  360 

8 

,0  iO 

3  x  4 

-,0  n  0 

7  x9 

26°  x  30° 

60°  x 100° | 

1 803  < 360  ' 

9 

3°  x 4° 

16°xi8° 

38°  x42° 

1 

60° x 100° ( 

1 

180°  x 360“ 

TABLE  2.1  Data  region  sizes  (A$>,A\),  centered  at  the 
computation  point. 


rms  (mgal) 

.  . 

5 '  x5  1 

15 ' x 15 ' 

-  I 

O 

»~4 

X 

O 

f— 4 

5°x5°  (fl4)  |  5  0  x  5  3  (DO 

a 

±8 

±8 

±5 

! 

I 

+  3  ;  -  3 

i 

b 

±10 

±7 

i 

±4 

i 

±3  ,  ±5 

! 

TABLE  2.2  Estimated  rms  errors  of  mean  anomalies. 


Error  correlations  were  to  be  neglected  throughout  (data  error  co- 
variance  matrix  is  diagonal) ,  a  questionable  assumption  in  the 
author's  opinion. 

Five  levels  of  prediction  (=  height  of  the  computation 
point)  were  considered: 
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hp  =  (30  000,  40  000,  70  000 ,  100  000,  200  000)  feet. 

The  estimated  quantity  is  the  gravity  disturbance  vector  5gp 
which  is  defined  as  the  difference  vector  between  the  actual  and 
the  normal  gravity  vector, 

5gp  =  gp  -  yp  ,  (2.i) 

taken  at  the  same  point  (Heiskanen  and  Moritz,  1967,  p.  85) .  The 
gravity  vector  is  the  gradient  of  the  corresponding  potential; 
therefore,  the  gravity  disturbance  vector  is  the  gradient  of  the 
disturbing  potential  T  , 


Sgp  =  gradpT 


(2.1)  • 


The  following  chapters  deal  with  accuracy  estimation  processes 
for  the  3  components  of  <5g  ,  the  radial  component  6  ,  and 

the  two  horizontal  components  (along  the  meridian)  and  3, 

(along  the  parallel) . 
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3 .  COLLOCATION  SOLUTION 


The  mean  square  error  of  a  predicted  gravity  field  quantity 
can  be  estimated  by 


m2  =  C 
p  pp 


-  C^C_1C 
p  p 


(3.1) 


(Moritz,  1980,  p.80) .  Here  m2  denotes  the  predicted  mean  square 
error,  Cpp  the  variance  of  the  predicted  quantity,  and  the  po¬ 
sitive  (strictly  speaking:  non-negative)  quantity  CpC-1Cp  the 
gain  of  prediction;  Cp  is  the  cross-covariance  vector  between 
the  predicted  quantity  and  all  data,  C  =  C  +  D  ,  where  C  denote 
the  data  covariance  matrix  and  D  the  a  priori  data  error  covaria 
ce  matrix.  D  reduces  to  diagonal  form  for  vanishing  error  corre¬ 
lations  . 

All  variances  and  covariances  are  derived  from  one  common 

covariance  function,  say  the  covariance  function  of  the  disturbing 

potential  K(P,Q)  .  K  is  harmonic  with  respect  to  P  and  Q 

outside  some  sphere  r  *  R_  ;  furthermore,  it  is  usually  chosen 

to  be  homogeneous  and  isotropic,  expressed  by  its  independence 

of  horizontal  position  and  direction.  This  is  why  K(P,Q)  depends 

only  on  the  product  of  the  moduli  of  the  radius  vectors  r  r 

R  y 

and  the  spherical  distance  between  the  two  points  P  and 

PQ 

Q  (which  are  located  outside  r  =  R  )  , 

B 


K  (P  ,Q)  « 


n  =  2 


I  n+  1 


r  r 

P  Q) 


pn(c°S'lp2) 


(3.2) 


(k  l  ,  n  =  2,  ...,  denote  (model)  decree  variances  and  P  the 
n  u 

n'th  degree  Legendre  polynomial. 

All  covariances  (which  enter  into  the  estimation  equation 
(3.1))  are  obtained  by  applying  the  rule  of  covariance  propatation 
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Since  the  linear  operator,  relating  the  disturbing  potential  to  'he 
gravity  anomaly ,  is  homogeneous  and  isotropic,  the  homogeneity  am 
isotropy  of  the  gravity  anomaly  covariance  function  follows  . 

A  least-squares  collocation  solution  is  burdened  by  two 
severe  problems:  the  calculation  of  the  many  individual  covariances 
and  the  inversion  cf  the  covariance  matrix.  Therefore,  every  pos¬ 
sibility  to  reduce  the  computational  effort  should  be  favoraoly 
considered.  In  the  cases  to  be  studied  here,  it  is  particularly 
the  regular  data  uistribution  which  car.  be  advantageously  taken 
into  account  in  the  structure  of  the  covariance  matrix  and  its  in¬ 
verse  . 

The  fact  that  mean  anomalies  are  to  be  used  as  data  requires 
special  consideration:  a  mean  anomaly  is  defined  over  3  "rectang¬ 
ular"  block  (i.',:')  ;  consequently,  the  linear  functionals  which 

have  to  be  applied  to  the  covariance  function  involves  an  integration 
over  that  particular  rectangle.  Because  cf  the  structure  of  the 
kernel,  (covariance  function)  ,  such  integrations  could  be  performed 
only  approximately  by  a  proper  numerical  integration  method,  a 
practically  impossible  enterprise:  first,  because  of  the  tremendous 
computational  effort,  and  second,  because  of  the  approximations 
involved.  (Integrating  the  covariance  function  numerically  over  a 
certain  region  means  approximating  the  covariance  function  by  a 
set  of  locally  restricted  polynomials,  usually  step  functions. 

Such  covariance  approximations  introduce  spectrum  disturbances  re¬ 
sulting  xu  partly  negative  eigenvalues  from  a  certain  degree  on 
(Siinkei,  19  78)  .  By  numerical  integration  of  the  covariance  function 
we,  therefore,  trade  in  an  integrated  effect  in  terms  of  possible 
singularities  cf  the  (anyway,  not  very  stable.)  covariance  matrix.) 

Fcr  these  reasons,  it  is  virtually  inevitable  to  replace  the  covari¬ 
ances  between  mean  values  at  zero  altitude  by  expressions  which, 
avoid  covariance  integrations.  Two  considerations,  which  go  back 
to  C.C.  Tscharning  (Tscherning  and  Rapp,  1974,  p.  70) ,  lead  to  the 
replacement  of  a  mean  anomaly  at  zero  altitude  by  a  corresponding 
point  anomaly  at  some  specified  altitude,  depending  on  the  block 
size  and  or,  the  parameters  of  the  covariance  function  :  a)  re- 


placing  the  rectangular  blocks  over  which  the  mean  anomalies  are 
defined  by  circular  regions  of  equal  si2e,  it  can  be  shown  that 
the  corresponding  mean  anomaly  covariance  function  is  obtained 
through  a  multiplication  of  the  degree  variances  by  the  square 
of  the  eigenvalues  of  the  moving  average  operator  (Scnwarz,  19',  6, 
p.  35  ff.) .  The  corresponding  infinite  series,  however,  can  no 
longer  be  represented  by  a  closed  covariance  expression;  the  in¬ 
finite  series,  could,  for  practical  purposes,  be  terminated  at 
a  certain  degree  n  -  200  for  5° *5'  ,  n  -  3000  for  5 '*5' 
block  size) ,  the  summations,  however,  are  too  time  consuming. 
Therefore,  there  seems  to  be  only  one  simple  way  to  overcome  this 
difficulty:  b)  replacing  the  squares  of  the  eigenvalues  of  the 
moving  average  operator  by  the  (n+1) ’ th  power  of  a  quantity 
u  <  1  ,  with  u  optimally  fitted.  un^ 1  can  be  easily  combined 
with  (R^/rpr^)n+1  which  leads  to  two  interpretations:  either  a 
diminishing  of  the  radius  of  the  Bjerhammar  sphere  R  or  an  up- 
ward  continuation  of  the  covariance  function  to  a  certain  altitude 
leaving  R  unchanged;  personally  I  find  the  latter  interpretation 
more  logical;  moreover,  it  preserves  consistency  in  the  domain  of 
homonicity.  The  u-values  and  altitudes  corresponding  to  different 
block  sizes  can  be  found  in  (Schwarz,  1976,  p.  39). 

Keeping  these  conceptual  replacements  in  mind,  we  can  from 
now  on  formally  consider  mean  gravity  anomalies  in  zero  altitude 
as  point  gravity  anomalies  in  a  certain  elevation.  In  the  following 
some  practical  and  theoretical  considerations  are  made  which  put 
emphasis  on  the  structure  of  the  covariance  matrix  for  gridded  data 
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3.1  Collocation  solution  for  gridded  data 


The  following  studies  deal  with  regularly  distributed  mean 
anomalies,  which,  for  reasons  explained  before,  can  be  treated 
like  gridded  point  values.  The  coverage  is  unfortunately  not 
total  (mean  anomalies  of  constant  block  sizes  are  restricted  to 
rectangular  zones  surrounding  the  computation  point)  .  A  total 
coverage  of  the  sphere  with  gridded  homogeneous  data  could  be 
treated  in  a  very  fast  and  elegant  way  taking  advantage  of  all 
existing  symmetries  (Colombo,  1979).  Even  in  the  case  of  a  partial 
gridded  data  coverage,  however,  the  structure  of  the  data  distribu¬ 
tion  carries  over  to  the  structure  of  the  corresponding  covariance 
matrix.  These  structures  can  advantageously  be  exploited  in  set¬ 
ting  up  and  inverting  the  covariance  matrix. 

Let  us  illustrate  such  a  typical  structure  by  means  of  a 
very  simple  example:  a  geographical  grid  of  data  with  4  rows 
(parallels)  and  6  colams  (meridians)  (Fig.  3.1).  Parallels  are 
assumed  to  be  equally  spaced  by  and  meridians  are  assumed  to 

be  equally  spaced  by  t 1  . 


FIG.  3.1  Gridded  data  arrangement 
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The  corresponding  covariance  matrix  is  a  6-4  x  6-4  matrix.  We 

partition  the  covariance  matrix  into  6x6  submatrices  each 

having  dimension  4x4  in  accordance  with  the  number  of  columns 

(6)  and  rows  (4)  .  The  arrangement  should  be  made  such  that  the 

column  counting  is  Mj ,  M2 ,  M3 ,  M_ x  ,  M_2,  M- ,  (from  inside  towards 

outside) .  If  we  denote  the  (4  x  4)  covariance  matrix  between 

two  columns  separated  by  mix  with  C  ,  the  full  covariance 

m 

matrix  has  the  following  structure: 


- 1 

n 

o 

Cl 

C2 

Cl 

1 

c  3 

Cl 

c0 

Cl 

C2 

C3 

c. 

C2 

C; 

Co 

C3 

C4 

C5 

Cl 

C2 

C3 

C  o 

1 

Cl 

C2 

C2 

C3 

c4 

Cl 

Co 

Cl 

C3 

- 

c4 

C5 

c2 

C, 

Co 

.J 

Three  properties  can  be  derived  from  (3.3)  immediately: 

a)  C  has  only  6  different  submatrices  Cg ,  . . . ,  C5  , 

b)  C  can  be  partitioned  into  2x2  blocks 


C  » 


E 

F 


F 

E 


with 


(3.4a) 


Co 

Cl 

c2 

Ci 

C2 

C31 

E  = 

C, 

Co 

Cl 

and 

F  = 

c2 

C3 

c4 ; 

C2 

1 

Cl 

c”. 

c3 

. 

C4 

C5i 

(3.4b) 


j 


) 

1 


c)  The  row-sum  of  C  equals  E  +  F  and  is  constant  -  a  very 
important  property  which  will  be  shown  to  account  for  sub¬ 
stantial  reductions  in  the  inversion  time. 

The  submatrix  Co  (covariance  matrix  of  column  M;  with  itself) 
is  structured  like  E  ,  the  other  submatrices  C.  ,  i  >  O 
(covariance  matrix  of  two  distinct  columns  separated  by  i  •  iU 
are  structured  like  F  ;  e.g. 


|  Cn  3 

f*0 

1 

2 

r •  0 

c0  3 

Coo 

Co  1 

Co  2 

Cc  3  ' 

1  „o 
|  C9  I 

^0  0 

c  ^ 

1 

c02 

Co 

Cl\ 

Cl  2 

1  1 

Cl  3  ! 

1 

|  cS2 

c  G 

1 

r*  ^ 
c0  0 

1 

/ 

Cl  = 

i 

L0  2 

Cl  2 

1 

C  2  2 

i  | 

C2  3  ! 

1  co 
'  c0  3 

Co  2 

^0  1 

c  o 

L0  0 

Co  3 

cj3 

1 

C2  3 

i  ! 
C  3  3  [ 

_ 

• 

L 

J 

(3.4c) 

Consequently,  we  observe  two  levels  of  structurization:  the  column- 
structure  (Co,  ...»  C 5 )  and  the  block-structure  (E,  F)  . 

Quadratic  matrices  with  constant  row-sum 

In  connection  with  probability  theory  we  will  sometimes 
find  matrices  for  which  the  sum  of  the  row-elements  is  constant 
and  equal  to  C  .  It  can  be  shown  (ZurmUhl,  1964,  p.  221  ff.) 
that  the  corresponding  inverse  has  also  a  constant  row-sum  equal 

to  1/C  . 

This  fact  leads  to  the  question  whether  a  generalization 
is  possible  such  that  the  elements  of  the  matrix  are  themselves 
submatrices . 
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Theorem:  Let  C  be  a  non-singular  quadratic  matrix  consisting 

of  quadratic  submatrices  of  equal  size,  and  let  the  row -sum  be 
constant  and  equal  to  E  .  Then  C  ^  consists  of  quadratic  sub¬ 
matrices  of  equal  size  and  has  a  constant  row-sum  equal  to  1  . 

Proof :  Obviously,  C  has  an  eigenvector-matrix  S  whose  elements 
are  unit-matrices,  and  E  is  the  eigenvalue-matrix, 


CS  =  S' 


(3.5a) 


It  follows  immediately  that  S  is  also  the  eigenvector-matrix  of 
C-1  ;  the  eigenvalue-matrix  is  E~ 1  : 


C  lS  =  SE 


(3.5b) 


Consequently,  E_l  is  the  row-sum  of  the  submatrices  of  C“ 1  .  o 


In  our  example  this  means  that 


f  E  F 


F  E 


G  H 


H  G 


E  =  E  +  F  , 


(3.6) 


G  +  H  =  (E  +  F) 


What  are  the  consequences  for  the  estimation  of  the  accuracy  of  a 
gravity  field  quantity  like  at  a  point  P  ,  located  on  the 

line  of  symmetry  with  respect  to  the  data  grid? 

Let  us  consider  the  cross-covariance  vector  Cp  ;  it  relates 
(statistically)  the  predicted  quantity  at  P  with  all  the  other 


data.  Since  the  data  are  regularly  distributed  on  a  grid,  whose 
line  of  symmetry  passes  the  prediction  point  P  ,  it  follows 
that  C  consists  of  two  subvectors  which  are  equal  to  each 
other. 


(C  is  the  cross-covariance  vector  between  the  predicted  quantity 
at  P  and  all  gridded  data  east  of  P  ;  it  is  at  the  same  time 
the  cross-covariance  vector  between  the  predicted  quantity  at  P 
and  all  gridded  data  west  of  P  .)  Introducing  (3.6)  and  (3.7) 
into  the  error  estimation  equation  (3.1),  we  obtain 


which  obviously  reduces  to 


m 


2  _ 


=  C. 


PP 


2  CP1<G 


+  H)  C 


PI 


or 


!m2 


=  C  -  2  CT, (E  +  F) _1C  . 
PP  PI  PI 


(3.8) 


(In  the  derivations  above  it  was  tacitly  assumed  that  the  data 
were  free  of  noise;  it  can  easily  be  shown,  however,  that  (3.8) 
holds  also  if  constant  noise  is  introduced  with  error  covariances 
dependent  on  the  spherical  distance.) 

Equation  (3.8)  shows,  how  one  single  symmetry  in  the  data 
pattern  (here:  symmetry  with  respect  to  the  meridian  of  the  pre- 
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r 


The  covariance  matrix  consists  of  4  submatrices  (2  in  the 
former  case) 


C1 1 

Cl2 

C  1  3 

C>4 

C  1  2 

Cn 

Cl4 

C1  3 

Ci3 

C  1  4 

C1 1 

C  1  2 

Cl  4 

C  1  3 

C12 

Cl  1 

(3.9) 


The  row-sum  is  obviously  constant  and  equal  to  l  with 

I  =  CU  +  C12  +  C13  +  C !  4  .  (3.10) 

In  analogy  to  (3.7)  the  cross-covariance  vector  between  the  pre¬ 
dicted  quantity  at  P  and  the  data  consists  of  4  equal  subvectors, 


Cp  = 


j  CP  l 
i  C 

L  P1J 


(3.11) 


Observing  the  above  theorem,  the  error  estimation  equation  (3.1) 
can  be  expressed  by 


m  =  C 
p  pp 


fcT  CT  CT  CT  ]  r  K 

[Si  'UP1  'Si  'Slj 


1 1 

K12 

K13 

K14 

12 

Ku 

K14 

*13 

1  3 

*14 

K1  1 

K12 

1  4 

*13 

K12 

Kn 

(3.12a) 


with 


Kn  +  K12  +  K,  3  +  K14  =  E-1  .  (3.12b) 

Consequently,  the  mean  square  error  is  simply  given  by 


m_ 


=  C 


PP 


-  4  C. 


:4C. 


+  c 


12 


+  C13  +  C 


-1, 


14 


'PI 


(3.13) 


It  should  again  be  noted  that  this  error  estimation  equation  is 
valid  if  the  data  pattern  has  two  lines  of  symmetry  with  respect 
to  the  prediction  point.  The  size  of  the  original  covariance 
matrix  to  be  inverted  is  split  into  4  ;  therefore,  a  gain  in 
speed  by  a  factor  of  64  can  be  expected  in  this  case  of  merid¬ 
ional  and  equatorial  symmetry. 


I. 


1 


INTEGRAL  APPROACH 


4  . 


It  is  well-known  that  the  disturbing  potential  T  can 
be  expressed  in  terms  of  free  air  anomalies  Ag  by 


T (r , $  ,X) 


S  ( r ,  41 )  A  gda 


with  S(r,<ii)  denoting  the  spatial  Stokes  kernel 


(4.1a) 


S  (r,!j/) 


2R  ,  R  ,  R1 
1  r 


r2 

COSl^ 


5  +  3  In 


r-Rcos'i'  +  l' 

2r 


(4.1b) 

(Heiskanen  and  Moritz,  1967,  p.  233). 

R  stands  for  the  mean  radius  of  the  earth,  \p  for  the  spherical 
distance  and  1  for  the  spatial  distance  between  two  points , 


1 


(r2  +  R2  -  2rRcos<|/) 


(4.1c) 


The  quantity  of  interest  is  the  gradient  of  T 


gradT  ==  |  |  gradpS (P ,Q) Ag (Q) da (Q)  , 


represented  in  terms  of  3  components 


«  (P)  «  £■ 

r  4ir 


D  S  (P  ,Q)  Ag  (Q)  da  (Q)  , 


J  r 


<5  .  (P) 


4^-ff  D*  S(P,Q)Ag(Q)da(Q)  , 


(4.2) 
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VP)  =  rffr  cfsV  Jj  °A  S(P,Q)Ag(Q)do(Q)  .  (4.2) 

p  p  1 J 

a  P 

If  the  gravity  anomalies  would  be  known  at  all  points  of 
a  and,  in  addition,  Ag  would  be  free  of  noise,  grad  T  could 
be  predicted  at  any  point  P  outside  the  surface  of  the  earth 
introducing  errors  only  due  to  spherical  approximation  and  due  to 
minor  conceptual  neglections  in  the  definition  of  Ag  ;  a 
detailed  discussion  can  be  found  in  (Heiskanen  and  Moritz,  1967, 
p.  240  ff .) . 

Naturally,  the  above  two  conditions  are  never  fulfilled: 
the  function  Ag  ,  postulated  by  (4.2),  is  usually  available  in 
terms  of  a  step- function  approximation  (=  mean  gravity  anomalies 
defined  on  rectangular  blocks);  in  addition,  these  mean  anomalies 
are  never  error-free.  Therefore,  two  additional  errors  are  intro¬ 
duced  in  the  calculation  of  grad  T  :  a  representation  error  (step 
function  representation  of  the  true  function  Ag) ,  and  a  data  error 
(upward  continued  integrated  effect  of  data-noise) .  The  following 
two  sections  deal  with  the  estimation  of  both  errors. 


4 . 1  Estimation  of  the  representation  error 

As  mentioned  before,  the  (unknown)  actual  function  Ag  is 
approximated  by  a  step  function;  the  size  of  the  steps  equal  the 
varying  size  of  the  blocks  (cf.  Table  2.1).  Representing  a  function 
(different  from  a  constant)  by  its  mean  values  means  loss  of  infor¬ 
mation,  particularly  in  the  higher  frequencies.  The  goal  is  to 
estimate  the  average  effect  of  such  information  deficiencies  onto 
the  gradient  of  the  disturbing  potential  in  high  altitudes.  Such 
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an  estimation  requires,  for  practical  reasons,  two  assumptions: 

a)  replacement  of  the  step  function  defined  over  rectangular 
blocks  by  a  moving  average  of  the  function,  such  that  the 
averaging  region  is  a  circle  with  an  area  equal  to  the  block 
area; 

b)  replacement  of  the  rectangular  zones  in  which  the  mean  values 
are  given  (cf.  Fig.  2.1)  by  circular  zones  of  equal  size. 

These  simplifications  are  introduced  in  order  to  obtain  isotropic 
operators;  note  that  similar  assumptions  have  been  made  for  the 

collocatioAa'iolution. 

. .«*•* 

A" 

Eigenvalues  of  integral  operators  with  isotropic  kernel 

Integral  operators  with  isotropic  kernel  play  a  fundamental 
role  in  representation  error  estimation  procedures . Therefore ,  a 
brief  sketch  of  basic  relations  will  be  given  in  the  sequel;  it 
is  essentially  an  outline  of  (Meissl,  1971b,  p.  38ff.). 

Consider  the  integral  transformation 

g  (P)  =  ||k(P,Q) f (Q)do (Q)  (4.3a) 

o 

with  an  isotropic  kernel  K  , 

K(P,Q)  =  K  ( cosijj )  .  (4.3b) 

.'.ccording  to  the  Funck-Hecke  formula  (Miiller,  1966  ,  p.20),  the 

spherical  harmonics  {$  }  are  eigenfunctions  of  this  integral 

nm 

transformation;  the  eigenvalues  <  are  projections  of  the  kernel 

n 

onto  the  Legendre  polynomials, 

f [k<P,Q)*  (Q)do(Q)  =  <  *n  (P) 

j  )  nm  n  nm 

a 


(4.4a) 


with 


l 

<  =  2it  [  K  ( t)  P  ( t)  dt  .  (4.4b) 

n  J  n 

-  1 

An  isotropic  integral  kernel  can  be  represented  in  terms  of  a 
series  of  Legendre  polynomials, 

K ( t)  =  I  k  P  (t)  ;  (4.5a) 

u  n  n 
n 

applying  (4.4b)  and  the  orthogonality  relation  of  Legendre  poly¬ 
nomials  , 


p  (t)p  (t)  = 
J  n  m 


it  follows  that  the  eigenvalues  of  K  are  given  by 


< 

n 


4tr 

2n  +  1 


k 


n 


(4.5b) 


Let 

f(P)  •  l  f  ♦  <P)  and  g  (P)  =  \  g  *  (P) 

nm  run  L  run  nm 

n  ,  m  n  ,  m 

be  spherical  harmonic  expansions  of  f  and  g  ,  respectively; 
then 


g  -  *  f  (4.6) 

nm  n  nm 

follows  from  (4.4a,b).  Equation  (4.6)  is  the  frequency  domain  ana¬ 
logue  of  the  integral  transformation  (4.3a).  It  will  play  a  central 
role  in  the  following. 


The  moving  average  operator 


If  the  isotropic  kernel  B(t)  of  an  integral  transformation 
like  (4.3a)  is  constant  for  t  i  t0  >  -1  and  vanishes  outside, 

B  is  called  a  moving  average  kernel, 


B  ( t ) 


1 

2tt  (l-t0  ) 


fOr  tg  i  t  Si  1 


0 


else 


(4.7) 


Its  eigenvalues  6  are  obtained  by  an  application  of  equation 

n 

(4.4b);  as  a  matter  of  fact,  they  depend  on  t0  , 


1 

S  (t0)  «  - - -  f  P  (t)  dt  . 

n  1  -  t0  J  n 

to 

The  following  expression  can  be  found  in  (Meissl,  1971a,  p.  24)  : 


Sn(tc 


1 


1  -  to  2n  +  1 


Pn+l(t0> 


(4.8) 


It  follows  from  a  Taylor-linearization  of  Pn(t)  at  t  =  1  that 

6  (1)  =  1  (...  function  reproduction).  (4.8)’ 

n 

Due  to  the  orthogonality  of  Legendre  polynomials  we  obtain 

B  (-1)  =  6  f  ...  function  anihilation)  . 
n  no 

Any  other  moving  average  operator  will  function  between  these  two 
extreme  cases,  reproduction  and  anihilation. 

It  is  now  a  simple  task  to  find  the  loss  of  information 
introduced  by  applying  a  moving  average  operator  onto  a  function  f , 


■***•.:  . 


5f  (P)  : 


f(P)  -  f(P)  =  l  [1 
n  ,  m 
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6  (tQ)  ]f  $  (P) 

n  u  nm  nm 


(4.9) 


The  spatial  Stokes  kernel  and  its  radial  derivative 

The  disturbing  potential  can  be  represented  by  a  series 
of  solid  spherical  harmonics 


T(r,0 , A) 


l 


n=«  2 


R 

r 


n+  1T 

n 


( 0  ,  A  )  ; 


(4.10) 


taking  into  account  the  relation  between  the  disturbing  potential 
T  and  the  gravity  anomaly  &g  in  the  frequency  domain 


T 

nm 


(4.11) 


T  can  be  formulated  in  terms  of  &g  , 


T(r, e , a)  =  R  l 


n+  1 


n 


- r  A9  (9  / 

l  n 


which,  by  expressing  the  Laplace  spherical  harmonic  a  9  (e  ,x  )  ex- 

n 

plicitly,  goes  over  into 


T (r ,  9,  A) 


_R_ 

4  TT 


V  R'n+I2n+1 

^  r  n- 1 
n=2  ' 


P^COStJ,) 


Agda 


The  infinite  sum  is  known  as  the  spatial  Stokes  function  S(r,^) 
which  is  given  in  its  closed  form  by  equation  (4.1b)  , 


S(r,*)  =  l  (|)n+1~r  .  (4.12) 

n  =  2V  1 

It  is  the  upward  continued  Stokes  kernel  S(R,ij>)  .  Its  radial 
derivative  follows  immediately. 


D^S  (r,t|/) 


n+1  (2n+l)  (n+1) 
n  -  1 


Pn  (cosi/<) . 


(4.13) 


— m«— i  t,  » 


yc •JWT*  »»—  -  -  ■ 
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Eigenvalues  of  the  kernel  - 


D  S  (r,'>) 


truncated  at 


In  order  to  estimate  the  impact  of  neglecting  Ag  -  infor¬ 
mation  in  circular  zones  around  the  computation  point  on  6r  ,  it 
is  essential  to  know  the  eigenvalue  of  the  corresponding  truncated 
kernel . 

If  the  kernel  -r^irD^S  ( r ,  ^ )  vanishes  for  '}><vo  >  its 
eigenvalues  are  obtained  by  applying  equation  (4.4b)  , 


sA(r't0) 


l 

k  =  2 


k+1 


( 2k+l )  (k+1) 
k  -  1 


Pn(t)Pk(t)dt.  (4, 


14) 


According  to  (Paul,  1973),  the  integral  of  products  of 
Legendre  functions  is  given  by 


R 


nk 


t(P 

. 

-  1 


P  (t) P  (t) dt 
n  k 


1 

(n-k) (n+k+1) 


n ( n+ 1 ) 


2n+l 


P.  (p  , 
k  n+ 1 


P 

n- 


1 


) 


k (k+1) 
2k+l 


(4.15a) 


for  n  ?  k  ,  and  the  recursion  formula 


R 


nn 


(n+1 )  ( 2n-l )  R  _  tH  2n-l  R 

n ( 2n+l )  n+l,n-l  n  n,n-2  2n+l  n-l,n-l 


(4.15b) 

for  n  »  k  ;  (here  and  in  the  sequel  we  suppress  the  argument  t0 
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for  the  sake  of  simplicity.)  Initial  values  are 


Roo  =  l  +  to 

Rn  =  (1  +  tJ)/3  .  (4.15c) 


What  follows,  is  a  tour  de  force,  the  attempt  to  find  a  closed 
expression  for  the  infinite  sum  in  equation  (4.14). 

Denoting  R/r  by  p  and  using  the  relations  (4.l5a,b), 
the  eigenvalues  s^  are  given  by 

3n(r't»)  =  j, 
k  =  2 

k^n 

Explicitly  written 


(2k+l) (k+1)  k+l„  .  ( 2n+l ) (n+1 )  n+l„ 

.  p  K  *r  “  ,  p  K 

k  -  1  nJc  n  -  1  nn 


(4.16a) 


s 


I 

n 


n(n+l)  . 

2n+l  'n+1 


P  i  V  ( 2k+l )  (k-H)  k+  lp 

n-r  £  (k-l)  (n-k)  (n+k+l)  p  k 

n>tk 


pP  l 
n  L 


k  =  1 
k/n-  1 


(k+2) 2 (k+1) 
k(n-k-l) (n+k+2) 


pk+1P, 


ip  r  _ k2  (k-1) _  k+1  ( 2n+l )  (n+1)  n+l 

p  n*  (k-2) (n-k+1) (n+k)  p  *k  n  1  p 

K  *  J 

k>*n  +  1 


(4.16b) 


We  know  (Tscherning,  19  72)  that  closed  expressions  are  available 
for  series  of  the  following  general  form: 
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F 

a 


(p  ,t)  : 


w 

y  _L_ 

+■  k+a 


k  =»0 


pk+1P,  (t) 
k 


and 


F 

a 


(p  ,  t) 


l 


k=  l-ct 


1 

k+a 


k+1 


P 


P  (t)  , 
k 


a  >  O 


a  £  0 


(4.17a) 


(4.17b) 


(Moritz,  1980,  p.  186)  where  a  denotes  a  fixed  integer.  There¬ 
fore,  we  express  the  coefficients  of  the  Legendre  functions  in 
terms  of  partial  fractions  of  the  form  (k+a)”1  ;  the  following 
identities  can  immediately  be  verified  : 


( 2k+l )  (k+1)  6 _ 1 _ n+1  1 

(k-1) (n-k) (n+k+l)  =  (n+£) (n-1)  k- 1  “  n-1  k-n 


n  1 
S+7  Ic+ri+T  ' 


(4.18a) 


(k+2)  2  (k+1 )  ^  4  1  n  (n+1)  2 _ 1 

k(n-k-l) (n+k+2)  (n+2) (n-1)  k  (2n+l)(n-l)  k-(n-l) 


+ 


n  2 ( n+ 1 )  1 

( 2 n+1 ) ( n+2)  k+n+2 


(4.18b) 


k2(k-l)  4 _ 1_  n(n+l)  2  1 

(k- .)  (n-k+1)  (n+k)  "  (n+2)  (n-1)  k-2  ~  (2n+l) (n-1)  k-(n+l) 


.  n  2  ( n+ 1 )  1 

( 2n+ 1 ) ( n+ 2 )  k+n 


1 


(4.18c) 


Introducing  these  partial  fractions  into  (4.16b)  leads  to  infinite 
sums  of  the  following  kind: 


h 

B 

c 

1 

r  1  k-t-1 

l  jT  »  Pk 

k=l  K  K 

y  -L.  ak+ip 
ki3  k'2  k 

k/n 

k  ^n-  1 

kj*n  +  1 

2 

y  1  0k+ip 

y  1  Jk+1r 

y  1  Dk+ir 

k«2k“n  k 

ki,  k-(n-l)  p  *k 

k- (n+1 )  p  Fk 

k^n 

k^n-  1 

k^n+  1 

3 

y  1  Dk+1r 

’  i  k+ir 

y  1  k+lp 

J2k+n+l  k 

ki1  k+n+2  p  k 

k“3  k+n  k 

k^n 

k^n-l 

k^n+  1 

I 

00 

p  k+ 1_ 

J,  0  * 

k*  1 

l  pk+1?k 

k=  3 

■ 

kj*n-  1 

_ 

kj*n+  1 

TABLE  4.1  Infinite  series  with  closed  expressions. 


Using  the  notation  of  (4.17a,b),  closed  expressions  are  available 
for  a  Z  -2  in  (Moritz,  1980,  p.  187,  188);  only  the  elements 
of  the  2nd  row  in  Table  4.1  pose  problems: 

Paul  (1973,  p.  418  ff.)  defines  a  quantity 
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U  (t,p)  can  also  be  related  to  F  defined  in  (4.17);  recalling 
n  a 

the  definition  (4.19)  ,  we  obtain 


F 

-  n 


n-  1 

I 


k-0 


1 

k-n 


and  after  straightforward  operations  involving  the  substitution 
of  (4.23),  we  finally  find  a  recurrence  relation  for  F  , 


II  —  ,  - 

I  s=r W 


n-1  .  1  k+ Ip  F  I 

n  [  ki0  k- (n-2)  p  Pk  F- ( n- 2 ) j 


-  A  "j1  ji-  „k*y  *  &  +  s^lsrj—IjL 

{  do  k“n  k  nj  n(2n-l) 


(4.24) 


With  the  initial  values 

Uq  ( t , p )  *  Fx  (t,p)  *  m  [1  +  , 

1  2' 

Uj (t,p)  =  7  F0 (t,p)  =  In  § 

p  t 

(N  :  *  1  +  L  -  pt)  , 


we  find  (Jn  for  arbitrary  n  by  applying  the  recursion  formula 
(4.23)  . 

Fn  with  n  >  0  and  a  corresponding  recurrence  relation 
can  be  found  in  (Moritz,  1980,  p.  188)  : 
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=  n  ( 1 )  (  p  -  P  )  I  _ — _  A  -  -  A  -  n  a 

(in+i)  n+1  n-l  [  (n+2)  (n-1)  1  n-1  2  n+2  A3_ 

+  ftP  f.  4 . a  .  n(n+l) 2  .  n2  (n+1 )  ' 

p  n  j_(n+2)  (n-1)  1  (2n+l)  (n-1)  a2  (2n+l)(n+2)  3  *4 

.  1  p  f  4  r  _  n(n+l) -  r  n2 (n+1) 

p  n  |_(n+2)  (n-1)  S  (2n+l)(n-l)  2  (2n+l)(n+2)  3 

r  1  .  (2n+l) (n+1)  n+ln 

4  n-1  p  nn 


(4.27) 


substituting  Ai#  Bi,  and  Ci  ,  the  eigenvalues  of  the  truncated 
kernel,  defined  in  (4.16a),  are  finally  given  by 
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(4.28) 


Again  the  dependence  of  to  has  been  suppressed  in  P,  U,  F,  R, 
and  s  . 

The  eigenvalues  s'  have  to  fulfil  various  conditions; 

n 

this  enables  us  to  roughly  check  the  equation: 

a)  s'(p,-l)  —  0  because  of  coinciding  limits  of  integration; 

.  .  _ . ,  , ,  _fR')n+1n+l  -  _  _  _  -  .. 


b)  s'(p,l)=  2  —t  — r-  for  n  2  2  ;  this  follows  from  the 

n  ( r  j  n- 1 

orthogonality  of  the  Legendre  polynomials; 


The  interested  reader  will  note  the  relation  of  s’  to  the  so¬ 
il 

called  Molodensky  coefficients  Q  ,  defined  as  the  eigenvalues 
of  the  isotropic  kernel  S{<|i)/2ir  ,  where  S('j/)  denotes  Stokes' 
function  (Heiskanen  and  Moritz,  1967,  p.  259  ff.), 


Qn(t0)  : 


r 


S(t)Pn(t)dt  ; 


for  t0  =  1  we  obtain. 


Q  (1)  «  rfr  • 

n  n- 1 


therefore  s'(p,1)  is  related  to  Q  (1)  by 
n  n 


S’(p,l)  =  ^  ( n+ 1 )  Q 

n  r  n 


2(1+Qn) 


(4.29) 


Formula  (4.28)  can  easily  be  programmed;  a  very  efficient  routine 
has  been  developed  which  calculates  1000  coefficients  for  arbitrary 
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p  •«  1  and  -1  -j  to  s  1  within  0.9  seconds  CPU- time  on  a 

Univac  1100.  (This  value  agrees  with  an  estimate  obtained  by 
Paul  (1973,  p.  422)  for  the  computation  of  the  Molodensky  coeffi¬ 
cients  using  an  analogous  procedure.) 

Now  we  come  back  to  the  estimation  of  the  representation 
error.  Let  the  free  air  gravity  anomaly  f  be  expanded  into  a 
series  of  orthogonal  spherical  harmonics, 


f(P»>  -  l  ' 

n  ,  m 

then  we  know  from  (4.6)  that  the  gravity  disturbance  series  expan¬ 
sion  can  simply  be  obtained  through  a  multiplication  of  the  coeffi¬ 
cients  by  the  eigenvalue  which  relates  both  quantities ,  the  gravity 
disturbance  at  a  fixed  altitude  and  the  gravity  anomaly;  denoting 
the  gravity  disturbance  by  6 r  ,  it  follows  with  (4.2)  and  (4.13) 
that 


!r<P)  l 


n  ,  m 


S'<0,l)f  m 
n  nm 


♦  (P  ) 


nm 


is  its  spherical  harmonic  expansion  expressed  at  the  altitude  of 
P  in  terms  of  gravity  anomaly  coefficients  {f  }  .  If  we  neglect 

nm 

gravity  information  outside  <ji  =  around  the  computation  point 
P  ,  an  error  in  5  (P)  is  introduced  which  can  simply  be  obtained 
by  replacing  the  eigenvalue  s^(p,l)  by  the  eigenvalue  of  the  cor¬ 
responding  truncated  kernel  s^(p,t0)  , 


e6  (P) 
r 


I  l 


n  ,  m 


s* (p,t0)f  ♦  (P) 

n  nm  nm 


If  in  the  same  region  the  actual  gravity  anomaly  function  is  replaced 
by  its  moving  average,  the  coefficients  are  to  be  replaced  by  the 


eigenvalue  difference  of  the  corresponding  moving  average 
operator 


e.  (P)  =  -  |  l  s'(o,t0)[l  -  8  (t0  )  3  f  $  (P)  . 

5  2  L  n  nu  nm  nm 

r  n  ,  m 

The  principle  should  be  clear:  each  integral  transformation  can 
be  considered  as  a  two-dimensional  convolution  on  the  sphere; 
it  is  well-known  that  a  convolution  of  two  functions  in  the  space 
domain  is  equivalent  to  a  multiplication  of  the  corresponding 
spectral  values  in  the  frequency  domain;  in  the  case  discussed 
above,  two  convolutions  are  performed  successively. 

In  the  case  to  be  studied  in  this  report,  5  regions 
are  defined  with  various  moving  averages.  Denoting  the  cosine  of 
the  radius  of  the  moving  average  circle  for  region  i  by  t  ^  and 
the  cosine  of  the  outer  radius  of  the  circular  zone,  in  which  this 
moving  average  is  applied,  by  t  ,  the  representation  error  is 
obtained  by  adding  the  corresponding  eigenvalue  products , 

e$  (P) 

r 

with 


=  -  %  y  y  s'  (P,t.  ,)[6(t.  ,)  -  s (x . ) ] f  *  (P) 

2  -  L  .  n  l-l  l-l  i  nm  nm 

n  ,  m  i=  1 

(4.30) 


0  (t)  =  6  (1)  =1  and  tc  =  1 
non 


I  stands  for  the  number  of  regions.  With 


Xn  :  =  2  j-  Sn(p,ti-l)(e(Ti-l)  '  S(Ti} 
l  =  1 


(4.31) 


we  can  estimate  the  mean  square  representation  error, 
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=  M{e2  (P)\  =  Mi  t  \'f  t  (P)  l  >' 
6  16  /  “  n  nm  rim  ^  r 

r  r  ;  v  n  ,  in  r  ,  s 


'f  $  (P) 

r  r  s  r  s 


which ,  due  to  the  orthogonality  of  the  basis  { $>  }  ,  reduces 

nm 


m2  =  )  \  - f 2 

:■  “  n  n  m 

r  n  ,  m 


!4.32) 


Recall  that  f  are  the  harmonic  coefficients  of  the  gravity 
nm 

anomaly's  expansion  into  a  series  of  fully  normalized  spherical 

harmonics.  Therefore,  the  sum  of  f2  taken  over  m  represents 

nm 

the  actual  gravity  anomaly  degree  variance  of  degree  n  , 

(Heiskanen  and  Moritz,  1967,  p.  259), 


c  :  =  l  f2 
n  L  nm 

m 


(4.33) 


Empirical  degree  variances  are  available  up  to  relatively  small  n 
(say  n  =  36)  ;  higher  degree  variances  have  to  be  taken  from  a 

degree  variance  model  like 


n  n+  2 

Ac0 


n  -  1 
(n-1) (n+B) 


n  >  2 


(4.34) 


with,  e.g. 


A  =  425  mgal2  , 

c  o  =  0.999617  , 

B  =  24  , 

c2  =  7.5  mgal2 

(Tscherning,  1976).  With  (4.33)  the  mean  square  representation 
error  is  expressed  by 


[ 

L 


oo 


(4.35) 


The  eigenvalues  sn(p,t0)  of  the  truncated  spatial  Stokes  kernel 
will  be  needed  in  order  to  derive  the  representation  error's 
impact  on  the  horizontal  components  of  the  gradient  of  the  disturbing 
potential.  With  (4.4b)  the  eigenvalues  of  the  kernel  S(r,*)/2* 
are  given  by 


s  (p,t0) 

n 


S ( r , t) P  (t) dt  , 
n 

-  1 


which,  by  substituting  equation  (4.12)  for  S(r,t),  is  expressed 
by 


sn (p ,t0 ) 


00 


-  l 

k  =  2 


2k+l 

k-1 


P  (t)Pu  ( t)  dt 
n  k 

-  1 


(4.36) 


The  integral  is  denoted  as  before;  it  can  be  expressed  in 

closed  form  in  terms  of  Legendre  polynomials  (equation  (4.15a) 
for  n  ^  k  ,  and  (4.15b)  for  n  =  k  )  .  With  this  notation  the 
infinite  series 


5  <P,t0)  *  l  .  + 

k=*2 
k/n 


,  •  ,  p  R 

nk  n-1  nn 


has  to  be  expressed  in  a  closed  form.  Introducing  (4.15a,b)  for 
and  representing  products  in  terms  of  partial  fractions 
(Paul,  1973,  p.  417,  418)  leads  to 
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(4.37) 


All  occurring  infinite  sums  have  already  been  used  before  and 
can  be  found  in  Table  4.1  ;  the  corresponding  closed  expressions 
are  given  in  equations  (4 . 26a ,b,c)  .  Introducing  them  in  (4.37), 
leads  to  the  explicit  form  of  the  eigenvalue. 
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(4.38) 
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The  eigenvalues  have  to  fulfil  the  following  conditions: 

a)  s  (p,-l)  =  0  because  of  coinciding  limits  of  integration, 

b)  sn(p  ,1)  =  pn+1  for  n  2  2  ;  this  follows  from  the 

orthogonality  of  the  Legendre  polynomials . 

As  a  matter  of  fact  s  ( 1 , to)  equals  the  Molodensky  coefficient  Q 

n  n 

defined  before, 

s  (l,tJ  =  Q  (tn)  . 
n  v  n  0 

Let  the  gravity  anomalies  f  be  expanded  into  a  series  of 
Laplace  surface  harmonics, 


f(P0)  =  l  fn(P0)  ; 
n 

then  we  know  from  (4.6)  that  the  disturbing  potential  at  the  level 
of  P  can  simply  be  obtained  through  a  multiplication  of  the  La- 
Dlace  surface  harmonics  by  the  eigenvalues  which  relate  both  quanti¬ 
ties;  it  follows  from  (4.1a)  that 

T (P)  =  |  l  sn(p ,l)f  (P)  .  (4.39) 

n 

If  gravity  information  is  neglected  outside  a  cap  of  radius  i|/  =  $0 
around  the  computation  point  P  ,  an  error  is  committed  which  can 
be  calculated  by  simply  replacing  the  eigenvalue  s  (p,1)  by  the 
eigenvalue  of  the  corresponding  truncated  kernel  s  (p,to)  , 


ST(P)  =  f  l  Sn(o't0)fn(P)  * 
n 

Following  identical  arguments  as  in  the  case  of  the  radial  derivative 
of  T  ,  and  using  the  same  notation  3  for  the  eigenvalues  of  the 
moving  average  kernel,  we  obtain  a  total  representation  error  in 
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the  disturbing  potential. 


eT(P)  =  f  l  l  -  S(r.)jf^P)  (4.40) 

n  i  =  1 

with  8  (t0)  =  6(1)  al  and  t0  =  1  ;  I  stands  for  the  number 
of  regions.  With 


Xn  :  =  I.^1Sn(p'ti-l)C6(Ti-l)  '  6(Ti} 
1=1 


(4.41) 


equation  (4.40)  reduces  to 


e  (P)  =  r  I  A  f  (P)  . 

T  n  n 

n 


(4.42; 


At  this  point  we  turn  from  the  disturbing  potential  to  the  horizontal 

components  of  its  gradient,  —  T  and  - — - -  D  T  ;the  gravity 

rp  $p  rpcos^p  Ap 

anomaly  representation  error  has  an  impact  onto  these  two  quantities 
which  can  directly  be  derived  from  (4.42)  : 


e6  (p) 


l  A  D*f 

L  n  $  n 
n 


(4.43) 


e.  (P) 

6a 


— —  y  a  d  .  f 

cos  ip  -  n  A  n 
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The  mean  square  error  of  the  total  horizontal  component 


me  =  M  H  +  es. 


is  thus  given  by 


f  ,  1 
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which,  according  to  (Heiskanen  and  Moritz,  1967,  p.  262),  reduces 
to 


i 


*  1  X2n(n+l)c 

e  %  n  n 

n  =  2 


(4.44) 


with  the  degree  variances  defined  in  (4.33). 

Equations  (4.35)  and  (4.44)  are  the  desired  representation 
error  estimates  for  the  radial  and  the  total  horizontal  component 
of  the  disturbing  potential's  gradient. 


4.2  Estimation  of  the  error  due  to  data  inaccuracy. 


In  the  foregoing  section  the  error  has  been  estimated  which 
resulted  from  a  replacement  of  the  (unknown)  actual  gravity  anomaly 
function  by  mean  values  with  variable  block  size  depending  on  its 
distance  from  the  computation  point  (see  Table  2.1);  in  order  to 
obtain  manageable  expressions,  the  mean  value  representation  was 
formally  replaced  by  a  corresponding  moving  average.  The  estimated 
error  has  been  called  representation  error. 

In  general,  the  mean  gravity  anomalies  are  affected  by  noise 
which  primarily  stems  from  the  mean  value  estimation  process.  This 
noise  is  the  second  important  error  source.  Its  impact  on  the 
gradient  of  the  disturbing  potential  could  be  calculated  if  the 
error  covariance  function  would  be  known;  this,  however,  is  hardly 


ever  the  case.  Under  such  circumstances,  it  is  a  widely  applied 
practice  to  consider  the  anomaly  estimates  as  independent. 

Performing  the  horizontal  differentiations  in  (4.2),  we 
obtain  the  gradient  of  T  in  the  following  form  (Heiskanen  and 
Moritz,  1967,  pp.  233,  234) 


Vp)  -  <7 
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^  PQ  U 


According  to  the  individual  block  sizes  of  the  mean  values,  the 
integral  has  to  be  splitted  up  into  J  subintegrals  if  J  mean 
anomalies  are  to  be  considered  for  the  estimation  of  the  grad  T  . 
Denoting  the  ms  estimate  of  the  j'th  mean  anomaly's  error  by  nu  , 
its  propagation  into  the  components  of  grad  T  is  as  follows 
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The  integrals  can  be  evaluated  using  a  fast  numerical  integration 
procedure  described  in  (SUnkel  and  Rummel,  1981). 
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5.  COLLOCATION  VERSUS  INTEGRAL  APPROACH 


Two  methods  of  gravity  disturbance  vector  estimation  have 
been  presented  here,  the  least-squares  collocation  solution  and 
the  integral  solution.  Although  they  seem  to  be  totally  different, 
they  share  a  number  of  common  features . 

Let  us  begin  with  the  data: 

The  accuracy  of  the  three  components  of  the  gravity  disturbance 
vector  is  estimated  on  the  basis  of  mean  free  air  gravity  anomalies; 
mean  anomalies  are  defined  on  "rectangular"  blocks  bounded  by  meri¬ 
dians  and  parallels,  therefore,  they  are  bound  to  the  global  coordi¬ 
nate  system;  the  operator  which  can  be  thought  of  transforming  the 
actual  gravity  anomaly  function  into  its  mean  value  representation 
is  non-isotropic. 

Covariance  functions,  which  are  commonly  in  use,  are  iso¬ 
tropic;  kernels  of  integral  formulas  are  also  mainly  isotropic.  An¬ 
isotropy  causes  problems  with  the  integration  of  covariance  functions 
and  with  the  integration  of  integral  kernels.  In  order  to  avoid 
these  difficulties,  isotropy  is  artificially  produced  through  a 
replacement  of  the  mean  value  concept  by  the  concept  of  moving 
average  over  a  circular  region.  Both  collocation  and  integral 
formulae  (representation  error  estimation)  use  the  eigenvalues  6 

n 

of  the  moving  average  kernel;  in  the  estimation  of  the  representa¬ 
tion  error,  the  eigenvalues  are  used  explicitly,  in  least-squares 
collocation,  they  enter  implicitly  via  the  covariance  function. 
(Usually,  an  approximation  is  used  in  order  to  obtain  closed  co- 
variance  expressions . ) 

The  data  error  estimates  enter  in  collocation  through  the 
error  covariance  matrix,  which  has  diagonal  form  if  no  correlations 
are  assumed;  the  same  diagonal  form  has  been  assumed  for  the  esti¬ 
mation  of  the  error  due  to  data  inaccuracy. 
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As  far  as  the  use  of  statistical  gravity  field  information 
is  concerned,  there  is  also  no  difference  between  the  two  methods: 
both  use  this  information  in  terms  of  a  degree  variance  model 
(and  empirical  degree  variances  for  lower  degrees) ;  again,  it 
enters  into  collocation  implicitly  via  the  covariance  function 
and  appears  explicitly  in  the  representation  error  estimation 
equation. 

The  upward  continuation  is  contained  in  the  eigenvalues 
of  the  kernel  (representation  error)  and  implicitly  in  the  indi¬ 
vidual  covariances  (collocation) . 

If  no  data  are  available,  both  methods  provide  identical 
error  estimates;  if  the  gravity  anomaly  would  be  known  at  every 
point  of  the  sphere,  the  collocation  solution  would  again  coincide 
with  the  integral  solution.  For  these  reasons,  we  can  expect  that 
the  estimations  will  differ  only  little  if  the  gravity  coverage  is 
reasonably  good  -  a  fact  which  was  strongly  confirmed  by  numerical 
calculations . 
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6.  RESULTS  -  CONCLUSIONS 


The  error  estimations  of  the  radial  and  horizontal  compo¬ 
nents  which  are  presented  here,  have  been  carried  out  by  least- 
squares  collocation  and  /  or  integral  formula  approach.  The  devia¬ 
tion  was  found  to  be  in  the  10%  range  and  less ,  dependent  main¬ 
ly  on  the  data  distribution;  therefore,  the  results  presented  here 
are  practically  valid  for  both  methods . 

The  numerical  investigation  was  started  with  a  small  data 
set  consisting  of  5'  x  5'  mean  anomalies  around  the  computation 
point;  between  4  and  144  anomalies  have  been  considered  for 
the  estimation;  two  cases  have  been  studied:  error-free  anomalies 
and  anomalies  with  a  rms-error  of  ±8  mgal;  the  prediction  point's 
latitude  was  assumed  to  be  zero;  Tscherning's  degree  variance  model 
2  has  been  used  with  36  lower  harmonics  subtracted;  this  corre¬ 
sponds  quite  well  to  the  assumption  of  error-free  5°  x  5°  anoma¬ 
lies  given  outside  the  small  region  in  which  5'  x  5'  anomalies 
are  available.  Of  course,  this  is  a  rather  poor  data  distribution 
and  is  not  representative  for  the  available  distribution;  I,  how¬ 
ever,  like  to  present  the  result  because  it  shows  the  essential 
behavior  of  all  solutions  remarkably  well: 

Radial  component: 

a)  The  error  decreases  rapidly  if  the  number  of  5'  x  5'  mean 
anomalies  around  the  computation  point  increases; 

b)  the  larger  the  5'  x  5'  data  set,  the  smaller  is  the  repre¬ 
sentation  error  (evident) ; 

c)  the  influence  of  data  errors  decreases  with  increasing  alti¬ 
tude; 

d)  a  very  typical  feature  is  the  cone-effect:  the  data  region 
of  strong  contribution  increases  with  increasing  prediction 
height,  or  with  other  words,  if  the  number  of  anomalies  is 


kept  constant,  the  prediction  error  increases  with  increasing 
altitude . 

Horizontal  component: 

a)  for  the  same  data  constellation,  the  error  is  more  than  twice 
as  high  and  decreases  only  relatively  slowly,  if  the  number 
of  5'  x  5'  anomalies  increases; 

b)  the  larger  the  data  set,  the  smaller  is  the  representation 
error  (evident) ; 

c)  the  influence  of  data  errors  decreases  with  increasing  alti¬ 
tude; 

d)  the  remote  zones  have  a  very  significant  influence  on  the  re¬ 
sult  (this  is  also  known  from  the  behavior  of  gravimetrically 
determined  deflections  of  the  vertical,  which  are  very  closely 
related  to  the  horizontal  components  of  grad  T) . 
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Table  6.1  contains  the  error  estimates  which  have  been  obtained 
using  the  data  set  described  in  Table  2.1  with  data  error  esti¬ 
mates  presented  in  Table  2.2.  From  these  figures,  we  can  draw 
the  following  conclusions: 

a)  For  a  fixed  altitude,  the  different  kinds  of  data  distribu¬ 
tion  (case  1  to  9  )  produce  no  significant  variation  in 
the  estimation  of  the  radial  component;  the  horizontal  compo¬ 
nent  estimates  differ  from  case  to  case  up  to  50%  . 

b)  The  rms-estimates  for  the  radial  component  are  strongly  in¬ 
fluenced  by  the  data  inaccuracy  which  accounts  for  up  to  90% 

of  the  total  error  in  low  altitudes;  the  data  error  impact  on  the 
the  radial  component  decreases  rapidly  with  increasing  eleva¬ 
tion.  Error-free  data  give  estimates  in  6^  of  ±0.4  mgal 
for  30  000  ft  down  to  ± 0.07  mgal  for  200  000  ft  altitude; 
these  figures  are  in  excellent  agreement  with  the  ones  derived 
from  the  representation  error  formula. 

With  the  available  data  (distribution,  rms-errors) ,  the 
radial  component  of  the  gravity  disturbance  vector  can  be  esti¬ 
mated  with  a  rms-error  of  ±1  mgal  at  an  altitude  of  about 
50  000  feet.  In  order  to  achieve  the  same  accuracy  in  30  000  ft 
elevation,  the  data  inaccuracy  has  to  be  reduced  by  about  60%  , 
particularly  in  the  region  around  the  computation  point  (5'  x  5' 
anomalies);  remote  zones  contribute  only  little. 

c)  The  situation  looks  much  worse  for  the  horizontal  components: 
compared  to  the  radial  component,  they  are  more  inaccurate  from 
a  factor  2  (at  30  000  feet  altitude)  up  to  a  factor  6  (at 
200  000  feet  altitude) ;  the  impact  of  data  inaccuracies  is 
relatively  little  (between  ±1.1  mgal  (30  000  feet)  and  ±0.2 
mgal  (200  000  feet));  the  crucial  point  is  the  data  distri¬ 
bution  ;  the  representation  error  varies  between  ±0.2  mgal 
(case  9)  and  ±2.9  mgal  (case  1)  and  decreases  very  slowly 
with  increasing  elevation  (less  than  10%  in  the  range  from 

30  000  to  200  000  feet) . 

With  the  available  data  (distribution,  rms-errors),  the 


48 


horizontal  components  of  the  gravity  disturbance  vector  can 
be  estimated  with  an  accuracy  of  ±2.3  mgal  at  30  000  ft 
elevation  with  the  best  data  distribution  available  (case  9a) ; 
this  corresponds  to  an  accuracy  of  ±0'.’5  for  the  direction 
of  the  gravity  vector.  An  improvement  to  ±1  mgal  (±0'.'2)  at 
30  OOO  ft  altitude  requires  a  considerably  better  represen¬ 
tation  of  the  gravity  anomaly  field;  unlike  in  the  case  of  the 
radial  component,  the  horizontal  component  responds  also  con¬ 
siderably  to  the  representation  in  the  medium  range  (up  to 
30°  spherical  distance  from  the  computation  point) .  In  this 
critical  region  the  block  sizes  need  to  be  reduced  by  a  factor 
of  about  2  and  the  overall  data  accuracy  should  be  increased 
by  about  30%  in  order  to  achieve  this  goal . 
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case 

30  000  ft 

radial  hor. 

40  000  ft 

radial  hor. 

70  000  ft 

radial  hor. 

100  000  ft 

radial  hor. 

200  000 

radial 

ft 

hor . 

la 

±1.5 

±3.1 

±1.2 

±3.0 

±0.7 

±2.9 

±0.6 

±2.9 

±0.4 

±2.8 

b 

±1.8 

±3.2 

±1.4 

±3.1 

±0.9 

±2.9 

±0.7 

±2.9 

±0.5 

±2.8 

2a 

±1.5 

±3.0 

±1.2 

±2.9 

±0.7 

±2.8 

±0.6 

±2.7 

±0.4 

±2.6 

b 

±1.8 

±3.1 

±1.4 

±2.9 

±0.9 

±2.8 

±0.7 

±2.7 

±0.5 

±2.6 

3a 

±1.5 

±2.8 

±1.2 

±2.7 

±0.7 

±2.6 

±0.6 

±2.6 

±0.4 

±2.5 

b 

±1.8 

±2.9 

±1.4 

±2.8 

±0.9 

±2.6 

±0.7 

±2.6 

±0.5 

±2.5 

4a 

±1.5 

±2.6 

±1.2 

±2.5 

±0.7 

±2.4 

±0.6 

±2.4 

±0.4 

±2.3 

b 

±1.8 

±2.8 

±1.4 

±2.6 

±0.9 

±2.4 

±0.7 

±2.4 

±0.5 

±2.3 

5a 

±1.5 

±2.8 

±1.2 

±2.7 

±0.7 

±2.6 

±0.6 

±2.6 

±0.4 

±2.5 

b 

±1.8 

±2.9 

±1.4 

±2.8 

±0.9 

±2.6 

±0.7 

±2.6 

±0.5 

6a 

±1.5 

±2.7 

±1.2 

±2.6 

±0.7 

±2.5 

±0.6 

±2.4 

±0.4 

±2.3 

b 

±1.8 

±2.8 

±1.4 

±2.6 

±0.9 

±2.5 

±0.7 

±2.4 

±0.5 

±2.3 

7a 

±1.5 

±2.7 

±1.2 

±2.6 

±0.7 

±2.5 

±0.6 

±2.5 

±0.4 

±2.4 

b 

±1.8 

±2.8 

±1.4 

±2.7 

±0.9 

±2.5 

±0.7 

±2.5 

±0.5 

±2.4 

8a 

±1.5 

±3.0 

±1.1 

±2.9 

±0.7 

±2.8 

±0.6 

±2.7 

±0.4 

±2.6 

b 

±1.8 

±3.1 

±1.4 

±2.9 

±0.9 

±2.8 

±0.7 

±2.7 

±0.4 

±2.6 

9a 

±1.5 

±2.3 

±1.1 

±2.2 

■H 

±0.6 

±2.0 

±0.4 

±1.9 

b 

±1.8 

±2.4 

±1.4 

±2.2 

■99 

B 

±0.6 

±2.0 

±0.4 

±1.9 

TABLE  6.1  rms-error  estimates  of  gravity  disturbance 
vector  components  in  various  altitudes 
(dimension:  mgal  ) 
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